---
title: "Bats and human blood project: Nyctalus Noctula torpor analysis"
output:
  html_document:
    df_print: paged
---

```{r setup, echo=FALSE}
knitr::opts_knit$set(root.dir = here::here())
# setwd(here::here())
```

# prepare & load packages
```{r}
pathWD <- "/Volumes/5360/Biophysics/Experiments/dRT-DC/220301_Bob_DataAnalysis_BatsHumans/20230823_NycNocTorporAnalysis"
setwd(pathWD)
# setwd(here::here())

pacman::p_load(tidyverse, lmerTest, car, emmeans, easystats, readr)

# load data ####
Bat1Table <- read_delim(file = "../Bat1Table_matProp.csv",
                        delim = ",",
                        col_names = TRUE)
Bat1Table$replicate <- paste0("NN", Bat1Table$index)

Bat2Table <- read_delim(file = "../Bat2Table_matProp.csv",
                        delim = ",",
                        col_names = TRUE)
Bat2Table$replicate <- paste0("NN", Bat2Table$index)

# concatenate data tables
combTable <- rbind(Bat1Table,Bat2Table)
remove(Bat1Table,Bat2Table)

# combTable
# View(combTable)
```


# define factors & sample size tables
```{r}
# define factors ####
combTable$animalType <- as.factor(combTable$animalType)
combTable$replicate <- as.factor(combTable$replicate)
combTable$tempLevel <- as.factor(combTable$tempLevel)
combTable$sex <- as.factor(combTable$sex)
combTable$operator <- as.factor(combTable$operator)
combTable$torporBlood <- as.factor(combTable$torporBlood)
combTable$torporCapture <- as.factor(combTable$torporCapture)

# levels ####
levels(combTable$animalType)
levels(combTable$replicate)
levels(combTable$tempLevel)
levels(combTable$sex)
levels(combTable$operator)

# overview table ####
writeLines("Torpor @ capture:")
table(combTable$tempLevel, combTable$torporCapture)
table(combTable$torporCapture)

writeLines("Torpor @ blood drawing:")
table(combTable$tempLevel, combTable$torporBlood)
table(combTable$torporBlood)

writeLines("Torpor @ capture:")
table(combTable$torporCapture, combTable$replicate)
writeLines("Torpor @ blood drawing:")
table(combTable$torporBlood, combTable$replicate)

```


# preparation
```{r}
levels.Temp <- levels(combTable$tempLevel)
levels.Torpor <- levels(combTable$torporBlood)

listOutcome <- c("area","tau1","tau2","deltaDefoTau1","deltaDefoTau2","eModulus","viscosity")
```

# density plots
```{r}
plots.Density <- purrr::map(listOutcome,
                           ~ ggplot(data = combTable,
                                    mapping = aes_(x = as.name(.x))) +
                             geom_density(size=1, alpha=0.3) +
                             ggtitle(paste(.x)))
plots.Density
```

# torporBlood
```{r}
# formulas
listFormula.tprBld.full <- purrr::map(listOutcome,
                                      ~ as.formula(paste(.x, "~ torporBlood * tempLevel + (1|replicate)"))
)

# calculate models
listModel.tprBld.full <- purrr::map(listFormula.tprBld.full,
                                    ~ lmer(formula = .x,
                                           data = combTable,
                                           REML = TRUE)
)

# step function from lmerTest package
# remove interactions which are not significant
# but keep always the main effects
listModel.tprBld.step <- purrr::map(listModel.tprBld.full,
                                    ~ lmerTest::step(.x,
                                                     ddf = "Satterthwaite",
                                                     alpha.fixed = 0.05,
                                                     alpha.random = 0.1,
                                                     reduce.fixed = TRUE,
                                                     reduce.random = FALSE,
                                                     keep = c("torporBlood", "tempLevel")
                             )
)

# extract final models
listModel.tprBld.final <- purrr::map(listModel.tprBld.step,
                                     ~ get_model(.x)
)

```


# torporCapture
```{r}
# formulas
listFormula.tprCptr.full <- purrr::map(listOutcome,
                                       ~ as.formula(paste(.x, "~ torporCapture * tempLevel + (1|replicate)"))
)

# calculate models
listModel.tprCptr.full <- purrr::map(listFormula.tprCptr.full,
                                     ~ lmer(formula = .x,
                                            data = combTable,
                                            REML = TRUE)
)

# step function from lmerTest package
# remove interactions which are not significant
# but keep always the main effects
listModel.tprCptr.step <- purrr::map(listModel.tprCptr.full,
                                     ~ lmerTest::step(.x,
                                                      ddf = "Satterthwaite",
                                                      alpha.fixed = 0.05,
                                                      alpha.random = 0.1,
                                                      reduce.fixed = TRUE,
                                                      reduce.random = FALSE,
                                                      keep = c("torporCapture", "tempLevel")
                             )
)

# extract final models
listModel.tprCptr.final <- purrr::map(listModel.tprCptr.step,
                                      ~ get_model(.x)
)

```


# print outome
```{r}
for (m in 1:length(listOutcome)) {
  writeLines("------------------------------------------------------------")
  writeLines(paste("Outcome variable:", listOutcome[m]))
  
  writeLines("-----Torpor @ blood drawing, Anova3-------------------------")
  # global p-values
  Anova(listModel.tprBld.final[[m]], type="3") %>% print()
  
  writeLines("-----Torpor @ capture, Anova3-------------------------------")
  # global p-values
  Anova(listModel.tprCptr.final[[m]], type="3") %>% print()
  
  writeLines("-----Torpor @ blood drawing, Summary------------------------")
  # summary (estimate) output
  summary(listModel.tprBld.final[[m]], corr=FALSE) %>% print()
  
  writeLines("-----Torpor @ capture, Summary------------------------------")
  # summary (estimate) output
  summary(listModel.tprCptr.final[[m]], corr=FALSE) %>% print()
  
  writeLines("-----End----------------------------------------------------\n")
  
}
```

# plot preparation
```{r}
library("ggh4x")
library("tools") # for capitalization

# labels for animal type
labels.Bool <- c("No", "Yes")
names(labels.Bool) <- levels.Torpor

# labels for temp
labels.Temp <- c("37", "23", "10")
names(labels.Temp) <- levels.Temp
labels.TempUnit <- paste(labels.Temp, "°C")
names(labels.TempUnit) <- levels.Temp

# axes titles
title.Torpor <- "Torpor"
title.Temp <- "Temperature"
title.TempUnit <- paste(title.Temp, "(°C)")
title.trpBld <- "Blood sampling"
title.trpCptr <- "Capturing"

# labels for outcome variables
listLabelOutcome <- c(bquote("Cell size" ~ ("μm"^2)),
                      expression(tau["inlet"] ~ "(ms)"),
                      expression(tau["channel"] ~ "(ms)"),
                      expression(hat(d)["inlet"]),
                      expression(hat(d)["channel"]),
                      "Young's modulus (kPa)",
                      "Viscosity (Pa s)")

colorPalette.Torpor <- c("#1F77BD","#8FBBDE")

```

# torpor, results
## torporBlood, results
```{r}
listResults.tprBld <- purrr::map(listModel.tprBld.final,
                                     ~ emmeans(.x,
                                               pairwise ~ torporBlood,
                                               lmer.df = "satterthwaite",
                                               adjust = NULL,
                                               lmerTest.limit = 1e6)
)

listResults.tprBld.means <- purrr::map(listResults.tprBld,
                                           ~ .x$emmeans %>% as.data.frame()
)
listResults.tprBld.means <- purrr::map(listResults.tprBld.means,
                                       ~ within(.x, {
                                         lower.SEM <- emmean - SE
                                         upper.SEM <- emmean + SE})
)

listResults.tprBld.contrasts <- purrr::map(listResults.tprBld,
                                               ~ .x$contrasts %>% as.data.frame()
)
```


## torporCapture, results
```{r}
listResults.tprCptr <- purrr::map(listModel.tprCptr.final,
                                     ~ emmeans(.x,
                                               pairwise ~ torporCapture,
                                               lmer.df = "satterthwaite",
                                               adjust = NULL,
                                               lmerTest.limit = 1e6)
)

listResults.tprCptr.means <- purrr::map(listResults.tprCptr,
                                           ~ .x$emmeans %>% as.data.frame()
)
listResults.tprCptr.means <- purrr::map(listResults.tprCptr.means,
                                       ~ within(.x, {
                                         lower.SEM <- emmean - SE
                                         upper.SEM <- emmean + SE})
)

listResults.tprCptr.contrasts <- purrr::map(listResults.tprCptr,
                                               ~ .x$contrasts %>% as.data.frame()
)
```


## plot preparation
```{r}
yMinFrac <- c(0.5, 0.2, 0.2, 0.5, 0.5, 0.7, 0.4)
sigFrac <- 0.2

# adapt contrasts data set for plotting
listResults.Torpor.means <- list()
listResults.Torpor.contrasts <- list()
listResults.Torpor.plotPara <- list()
for (m in 1:length(listModel.tprCptr.final)) {

  dataMax <- max( listResults.tprCptr.means[[m]] %>% select(upper.SEM) %>% max(., na.rm = TRUE), listResults.tprBld.means[[m]] %>% select(upper.SEM) %>% max(., na.rm = TRUE) )
  yLo <- dataMax * yMinFrac[m]
  yHi <- dataMax * (1 - sigFrac*yMinFrac[m]) / (1 - sigFrac)
  
  listResults.Torpor.plotPara[[m]] <- list(dataMax = dataMax,
                                               yLo = yLo,
                                               yHi = yHi)
  
  # means
  a <- listResults.tprBld.means[[m]]
  b <- listResults.tprCptr.means[[m]]
  a <- rename(a, torporFlag = torporBlood)
  b <- rename(b, torporFlag = torporCapture)
  
  a <- within(a, torporType <- title.trpBld)
  b <- within(b, torporType <- title.trpCptr)
  
  listResults.Torpor.means[[m]] <- rbind(a, b)
  
  # contrats
  a <- listResults.tprBld.contrasts[[m]]
  b <- listResults.tprCptr.contrasts[[m]]
  a <- within(a, torporType <- title.trpBld)
  b <- within(b, torporType <- title.trpCptr)
  
  listResults.Torpor.contrasts[[m]] <- rbind(a, b)
  
  
  listResults.Torpor.contrasts[[m]] <- listResults.Torpor.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value < 0.001 ~ "***",
                               p.value < 0.01  ~ "**",
                               p.value < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = c(1) + 0.05,
           xEnd     = c(2) - 0.05,
           yLineRel = (1 + c(0) ) / 1.8,
           yTextRel = 0.015,
           yStart = (yHi - dataMax)*yLineRel + dataMax,
           yEnd   = yStart,
           yText  = (yLineRel + yTextRel)*(yHi - dataMax) + dataMax,
           xText  = 0.5*(xStart+xEnd))

  rm(dataMax, yLo, yHi, a, b)

}
```


## all plots Torpor
```{r fig.height=8/2.54, fig.width=9/2.54}
plots.Torpor <- list()
tbl.plotResults.Torpor <- data.frame(matrix(ncol = 5, nrow = 0))

for (m in 1:length(listModel.tprBld.final)) {
  
  dataMax <- listResults.Torpor.plotPara[[m]]$dataMax
  yLo <- listResults.Torpor.plotPara[[m]]$yLo
  yHi <- listResults.Torpor.plotPara[[m]]$yHi
  
  plots.Torpor[[m]] <- ggplot(data = listResults.Torpor.means[[m]])
  
  # define layout parameters
  plots.Torpor[[m]] <- plots.Torpor[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          strip.text.y.right = element_text(size = 11, angle = 90),
          strip.text.x = element_text(size = 11),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none") +
  facet_grid(cols = vars(torporType))
  
  # create plot
  plots.Torpor[[m]] <- plots.Torpor[[m]] +
    geom_bar(mapping = aes(x = torporFlag,
                           y = emmean,
                           fill = torporFlag),
             stat = "identity", # important since we provide mean values already
             size = 0.75, # line width of bar borders
             colour = "black", # bar border color
             fill = rep(colorPalette.Torpor, times = 2),
             position = position_dodge()) # important for grouped bars
  
  # add error bars
  plots.Torpor[[m]] <- plots.Torpor[[m]] +
    geom_errorbar(mapping = aes(x = torporFlag,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  position = position_dodge())
  
  # add nested x-axis
  plots.Torpor[[m]] <- plots.Torpor[[m]] +
    scale_x_discrete(labels = labels.Bool)
  
  # add p-values
  plots.Torpor[[m]] <- plots.Torpor[[m]] +
    geom_segment(data = listResults.Torpor.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.Torpor.contrasts[[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  plots.Torpor[[m]] <- plots.Torpor[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  plots.Torpor[[m]] <- plots.Torpor[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, yHi),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(subtitle = title.Torpor,
         x = title.Torpor,
         y = listLabelOutcome[m],
         caption = paste0("Models: ", format(listFormula.tprBld.full[[m]]), "\n", format(listFormula.tprCptr.full[[m]])))
  
  # change coordinates = "zoom in"
  plots.Torpor[[m]] <- plots.Torpor[[m]] +
    coord_cartesian(ylim = c(yLo, yHi))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_Torpor", ".svg")),
         width = 9, height = 8, units = "cm", dpi = 600, device = "svg",
         plot = plots.Torpor[[m]], scale = 1)

  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_Torpor", ".png")),
         width = 9, height = 8, units = "cm", dpi = 600, device = "png",
         plot = plots.Torpor[[m]], scale = 1)
  
  # results table
  tbl.plotResults.Torpor <- rbind(tbl.plotResults.Torpor, data.frame(rep(listOutcome[m], times = nrow(listResults.Torpor.means[[m]])), listResults.Torpor.means[[m]]$torporType, listResults.Torpor.means[[m]]$torporFlag, listResults.Torpor.means[[m]]$emmean, listResults.Torpor.means[[m]]$SE, listResults.Torpor.means[[m]]$upper.CL - listResults.Torpor.means[[m]]$emmean))
}

colnames(tbl.plotResults.Torpor) <- c("targetVar", "torporType", "torporFlag ", "mean", "SEM", "CI")
write.csv(tbl.plotResults.Torpor, file.path(pathWD, paste0("tblPlotResults", "_Torpor", ".csv")), row.names = FALSE)

plots.Torpor

```


# torpor | temp, results
## torporBlood | temp, results
```{r}
listResults.tprBldTemp <- purrr::map(listModel.tprBld.final,
                                     ~ emmeans(.x,
                                               pairwise ~ torporBlood | tempLevel,
                                               lmer.df = "satterthwaite",
                                               adjust = NULL,
                                               lmerTest.limit = 1e6)
)

listResults.tprBldTemp.means <- purrr::map(listResults.tprBldTemp,
                                           ~ .x$emmeans %>% as.data.frame()
)
listResults.tprBldTemp.means <- purrr::map(listResults.tprBldTemp.means,
                                           ~ within(.x, {
                                             lower.SEM <- emmean - SE
                                             upper.SEM <- emmean + SE})
)

listResults.tprBldTemp.contrasts <- purrr::map(listResults.tprBldTemp,
                                               ~ .x$contrasts %>% as.data.frame()
)
```

## torporCapture, results
```{r}
listResults.tprCptrTemp <- purrr::map(listModel.tprCptr.final,
                                     ~ emmeans(.x,
                                               pairwise ~ torporCapture | tempLevel,
                                               lmer.df = "satterthwaite",
                                               adjust = NULL,
                                               lmerTest.limit = 1e6)
)

listResults.tprCptrTemp.means <- purrr::map(listResults.tprCptrTemp,
                                           ~ .x$emmeans %>% as.data.frame()
)
listResults.tprCptrTemp.means <- purrr::map(listResults.tprCptrTemp.means,
                                           ~ within(.x, {
                                             lower.SEM <- emmean - SE
                                             upper.SEM <- emmean + SE})
)

listResults.tprCptrTemp.contrasts <- purrr::map(listResults.tprCptrTemp,
                                               ~ .x$contrasts %>% as.data.frame()
)
```


## plot preparation
```{r}
yMinFrac <- c(0.5, 0.2, 0.2, 0.5, 0.5, 0.38, 0.3)
sigFrac <- 0.25

# adapt contrasts data set for plotting
listResults.TorporTemp.means <- list()
listResults.TorporTemp.contrasts <- list()
listResults.TorporTemp.plotPara <- list()
for (m in 1:length(listModel.tprCptr.final)) {

  a <- listResults.tprCptrTemp.means[[m]] %>% select(upper.SEM) %>% max(., na.rm = TRUE)
  b <- listResults.tprBldTemp.means[[m]] %>% select(upper.SEM) %>% max(., na.rm = TRUE)
  dataMax <- max(a,b)
  yLo <- dataMax * yMinFrac[m]
  yHi <- dataMax * (1 - sigFrac*yMinFrac[m]) / (1 - sigFrac)
  
  listResults.TorporTemp.plotPara[[m]] <- list(dataMax = dataMax,
                                               yLo = yLo,
                                               yHi = yHi)
  
  # means
  a <- listResults.tprBldTemp.means[[m]]
  b <- listResults.tprCptrTemp.means[[m]]
  a <- rename(a, torporFlag = torporBlood)
  b <- rename(b, torporFlag = torporCapture)
  
  a <- within(a, torporType <- title.trpBld)
  b <- within(b, torporType <- title.trpCptr)
  
  listResults.TorporTemp.means[[m]] <- rbind(a, b)
  
  # contrats
  a <- listResults.tprBldTemp.contrasts[[m]]
  b <- listResults.tprCptrTemp.contrasts[[m]]
  a <- within(a, torporType <- title.trpBld)
  b <- within(b, torporType <- title.trpCptr)
  
  listResults.TorporTemp.contrasts[[m]] <- rbind(a, b)
  
  
  listResults.TorporTemp.contrasts[[m]] <- listResults.TorporTemp.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value < 0.001 ~ "***",
                               p.value < 0.01  ~ "**",
                               p.value < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = c(1) + 0.05,
           xEnd     = c(2) - 0.05,
           yLineRel = (1 + c(0) ) / 1.8,
           yTextRel = 0.015,
           yStart = (yHi - dataMax)*yLineRel + dataMax,
           yEnd   = yStart,
           yText  = (yLineRel + yTextRel)*(yHi - dataMax) + dataMax,
           xText  = 0.5*(xStart+xEnd))

  rm(dataMax, yLo, yHi, a, b)

}
```


## all plots Torpor | torpor type, temp
```{r fig.height=14/2.54, fig.width=9.6/2.54}
plots.TorporTypeTemp <- list()
tbl.plotResults.TorporTypeTemp <- data.frame(matrix(ncol = 6, nrow = 0))

for (m in 1:length(listModel.tprBld.final)) {
  
  dataMax <- listResults.TorporTemp.plotPara[[m]]$dataMax
  yLo <- listResults.TorporTemp.plotPara[[m]]$yLo
  yHi <- listResults.TorporTemp.plotPara[[m]]$yHi
  
  plots.TorporTypeTemp[[m]] <- ggplot(data = listResults.TorporTemp.means[[m]])
  
  # define layout parameters
  plots.TorporTypeTemp[[m]] <- plots.TorporTypeTemp[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          strip.text.y.right = element_text(size = 11, angle = 90),
          strip.text.x = element_text(size = 11),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none")+
  facet_grid(cols = vars(torporType), rows = vars(tempLevel), labeller = labeller(tempLevel = labels.TempUnit))
  
  # create plot
  plots.TorporTypeTemp[[m]] <- plots.TorporTypeTemp[[m]] +
    geom_bar(mapping = aes(x = torporFlag,
                           y = emmean,
                           fill = torporFlag),
             stat = "identity", # important since we provide mean values already
             size = 0.75, # line width of bar borders
             colour = "black", # bar border color
             fill = rep(colorPalette.Torpor, 2*length(levels.Temp)),
             position = position_dodge()) # important for grouped bars
  
  # add error bars
  plots.TorporTypeTemp[[m]] <- plots.TorporTypeTemp[[m]] +
    geom_errorbar(mapping = aes(x = torporFlag,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  position = position_dodge())
  
  # add nested x-axis
  plots.TorporTypeTemp[[m]] <- plots.TorporTypeTemp[[m]] +
    scale_x_discrete(labels = labels.Bool)
  
  # add p-values
  plots.TorporTypeTemp[[m]] <- plots.TorporTypeTemp[[m]] +
    geom_segment(data = listResults.TorporTemp.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.TorporTemp.contrasts [[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  plots.TorporTypeTemp[[m]] <- plots.TorporTypeTemp[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  plots.TorporTypeTemp[[m]] <- plots.TorporTypeTemp[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, yHi),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(subtitle = title.Torpor, 
         x = title.Torpor,
         y = listLabelOutcome[m],
         caption = paste0("Models: ", format(listFormula.tprBld.full[[m]]), "\n", format(listFormula.tprCptr.full[[m]])))
  
  # change coordinates = "zoom in"
  plots.TorporTypeTemp[[m]] <- plots.TorporTypeTemp[[m]] +
    coord_cartesian(ylim = c(yLo, yHi))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_TorporTypeTemp", ".svg")),
         width = 9.6, height = 14, units = "cm", dpi = 600, device = "svg",
         plot = plots.TorporTypeTemp[[m]], scale = 1)
  
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_TorporTypeTemp", ".png")),
         width = 9.6, height = 14, units = "cm", dpi = 600, device = "png",
         plot = plots.TorporTypeTemp[[m]], scale = 1)
  
  # results table
  tbl.plotResults.TorporTypeTemp <- rbind(tbl.plotResults.TorporTypeTemp, data.frame(rep(listOutcome[m], times = nrow(listResults.TorporTemp.means[[1]])), listResults.TorporTemp.means[[m]]$torporType, listResults.TorporTemp.means[[m]]$torporFlag, listResults.TorporTemp.means[[m]]$tempLevel,  listResults.TorporTemp.means[[m]]$emmean,  listResults.TorporTemp.means[[m]]$SE, listResults.TorporTemp.means[[m]]$upper.CL - listResults.TorporTemp.means[[m]]$emmean))
}

colnames(tbl.plotResults.TorporTypeTemp) <- c("targetVar", "torporType", "torporFlag", "tempLevel", "mean", "SEM", "CI")
write.csv(tbl.plotResults.TorporTypeTemp, file.path(pathWD, paste0("tblPlotResults", "_TorporTypeTemp", ".csv")), row.names = FALSE)

plots.TorporTypeTemp

```


# temp | torpor, results
## temp | torporBlood, results
```{r}
listResults.TempTprBld <- purrr::map(listModel.tprBld.final,
                                     ~ emmeans(.x,
                                               pairwise ~ tempLevel | torporBlood,
                                               lmer.df = "satterthwaite",
                                               adjust = NULL,
                                               lmerTest.limit = 1e6)
)

listResults.TempTprBld.means <- purrr::map(listResults.TempTprBld,
                                           ~ .x$emmeans %>% as.data.frame()
)
listResults.TempTprBld.means <- purrr::map(listResults.TempTprBld.means,
                                           ~ within(.x, {
                                             lower.SEM <- emmean - SE
                                             upper.SEM <- emmean + SE})
)

listResults.TempTprBld.contrasts <- purrr::map(listResults.TempTprBld,
                                               ~ .x$contrasts %>% as.data.frame()
)
```

## temp | torporCapture, results
```{r}
listResults.TempTprCptr <- purrr::map(listModel.tprCptr.final,
                                     ~ emmeans(.x,
                                               pairwise ~ tempLevel | torporCapture,
                                               lmer.df = "satterthwaite",
                                               adjust = NULL,
                                               lmerTest.limit = 1e6)
)

listResults.TempTprCptr.means <- purrr::map(listResults.TempTprCptr,
                                           ~ .x$emmeans %>% as.data.frame()
)
listResults.TempTprCptr.means <- purrr::map(listResults.TempTprCptr.means,
                                           ~ within(.x, {
                                             lower.SEM <- emmean - SE
                                             upper.SEM <- emmean + SE})
)

listResults.TempTprCptr.contrasts <- purrr::map(listResults.TempTprCptr,
                                               ~ .x$contrasts %>% as.data.frame()
)
```

## plot preparation
```{r}
yMinFrac <- c(0.5, 0.2, 0.2, 0.5, 0.5, 0.2, 0.1)
sigFrac <- 0.35

# adapt contrasts data set for plotting
listResults.TempTorpor.means <- list()
listResults.TempTorpor.contrasts <- list()
listResults.TempTorpor.plotPara <- list()
for (m in 1:length(listModel.tprCptr.final)) {

  a <- listResults.TempTprCptr.means[[m]] %>% select(upper.SEM) %>% max(., na.rm = TRUE)
  b <- listResults.TempTprBld.means[[m]] %>% select(upper.SEM) %>% max(., na.rm = TRUE)
  dataMax <- max(a,b)
  yLo <- dataMax * yMinFrac[m]
  yHi <- dataMax * (1 - sigFrac*yMinFrac[m]) / (1 - sigFrac)
  
  listResults.TempTorpor.plotPara[[m]] <- list(dataMax = dataMax,
                                               yLo = yLo,
                                               yHi = yHi)
  
  # means
  a <- listResults.TempTprBld.means[[m]]
  b <- listResults.TempTprCptr.means[[m]]
  a <- rename(a, torporFlag = torporBlood)
  b <- rename(b, torporFlag = torporCapture)
  a <- within(a, torporType <- title.trpBld)
  b <- within(b, torporType <- title.trpCptr)
  
  listResults.TempTorpor.means[[m]] <- rbind(a, b)
  
  # contrats
  a <- listResults.TempTprBld.contrasts[[m]]
  b <- listResults.TempTprCptr.contrasts[[m]]
  a <- rename(a, torporFlag = torporBlood)
  b <- rename(b, torporFlag = torporCapture)
  a <- within(a, torporType <- title.trpBld)
  b <- within(b, torporType <- title.trpCptr)
  
  listResults.TempTorpor.contrasts[[m]] <- rbind(a, b)
  
  
  listResults.TempTorpor.contrasts[[m]] <- listResults.TempTorpor.contrasts[[m]] %>%
    mutate(sig.idx = case_when(p.value < 0.001 ~ "***",
                               p.value < 0.01  ~ "**",
                               p.value < 0.05  ~ "*",
                               TRUE ~ "")) %>%
    mutate(p.valuef = sprintf("%1.4f", p.value)) %>%
    mutate(p.valuef = case_when(p.valuef =="0.0000" ~ "< 0.0001",
                                TRUE ~ as.character(p.valuef))) %>%
    mutate(xStart   = rep(c(1,1,2) + 0.05, times = length(levels.Torpor)*length(levels.Torpor)),
           xEnd     = rep(c(2,3,3) - 0.05, times = length(levels.Torpor)*length(levels.Torpor)),
           yLineRel = rep((1 + c(0,1,0) ) / 2.8, times = length(levels.Torpor)*length(levels.Torpor)),
           yTextRel = 0.015,
           yStart = (yHi - dataMax)*yLineRel + dataMax,
           yEnd   = yStart,
           yText  = (yLineRel + yTextRel)*(yHi - dataMax) + dataMax,
           xText  = 0.5*(xStart+xEnd))

  rm(dataMax, yLo, yHi, a, b)

}
```


## all plots temp | torpor
```{r fig.height=11/2.54, fig.width=12/2.54}
plots.TempTorpor <- list()
tbl.plotResults.TempTorpor <- data.frame(matrix(ncol = 6, nrow = 0))

for (m in 1:length(listModel.tprBld.final)) {
  
  dataMax <- listResults.TempTorpor.plotPara[[m]]$dataMax
  yLo <- listResults.TempTorpor.plotPara[[m]]$yLo
  yHi <- listResults.TempTorpor.plotPara[[m]]$yHi
  
  plots.TempTorpor[[m]] <- ggplot(data = listResults.TempTorpor.means[[m]])
  
  # define layout parameters
  plots.TempTorpor[[m]] <- plots.TempTorpor[[m]] +
    theme_bw() +
    theme(axis.title.x = element_text(size = 15, margin = margin(t=10)),
          axis.text.x = element_text(size = 10, face = "bold", margin = margin(t=5)),
          axis.title.y = element_text(size = 15, face = "bold", margin = margin(r=10)),
          axis.text.y = element_text(size = 10, face = "bold"),
          strip.text.y.right = element_text(size = 11, angle = 90),
          strip.text.x = element_text(size = 11),
          panel.border = element_rect(size = 1.5, colour = "black"),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          plot.title = element_text(size = 20, face = "bold", hjust = 0.5, vjust = 1.0),
          plot.subtitle = element_text(size = 15, hjust = 0.5, vjust = 0),
          plot.caption = element_text(size = 8, face = "italic", vjust = 0, margin = margin(t=10)),
          legend.position = "none")+
  facet_grid(cols = vars(torporType), rows = vars(torporFlag), labeller = labeller(torporFlag = labels.Bool))
  
  # create plot
  plots.TempTorpor[[m]] <- plots.TempTorpor[[m]] +
    geom_bar(mapping = aes(x = tempLevel,
                           y = emmean,
                           fill = tempLevel),
             stat = "identity", # important since we provide mean values already
             size = 0.75, # line width of bar borders
             colour = "black", # bar border color
             position = position_dodge()) # important for grouped bars
  
  # add error bars
  plots.TempTorpor[[m]] <- plots.TempTorpor[[m]] +
    geom_errorbar(mapping = aes(x = tempLevel,
                                ymin = lower.SEM,
                                ymax = upper.SEM),
                  size=0.75,
                  width=0.25,
                  position = position_dodge())
  
  # add nested x-axis
  plots.TempTorpor[[m]] <- plots.TempTorpor[[m]] +
    scale_x_discrete(labels = labels.Temp)
  
  # add p-values
  plots.TempTorpor[[m]] <- plots.TempTorpor[[m]] +
    geom_segment(data = listResults.TempTorpor.contrasts[[m]],
                 mapping = aes(x = xStart, xend = xEnd, y = yStart, yend = yEnd), 
                 size = 0.5, colour="black", linetype="solid", inherit.aes = FALSE) +
    geom_text(data = listResults.TempTorpor.contrasts [[m]],
              mapping = aes(x = xText, y = yText, label = sig.idx),
              size = 7)
  
  # add nicer fill color
  plots.TempTorpor[[m]] <- plots.TempTorpor[[m]] +
    scale_fill_brewer(palette = "Blues")
  
  # add further layout parameter
  plots.TempTorpor[[m]] <- plots.TempTorpor[[m]] +
    scale_y_continuous(guide = "axis_minor", # guide type comes from package "ggh4x"
                       limits = c(0, yHi),
                       expand = expansion(mult = c(0.025, 0.025))) +
    labs(subtitle = title.Torpor, 
         x = title.TempUnit,
         y = listLabelOutcome[m],
         caption = paste0("Models: ", format(listFormula.tprBld.full[[m]]), "\n", format(listFormula.tprCptr.full[[m]])))
  
  # change coordinates = "zoom in"
  plots.TempTorpor[[m]] <- plots.TempTorpor[[m]] +
    coord_cartesian(ylim = c(yLo, yHi))
  
  # save plots
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_TempTorpor", ".svg")),
         width = 12, height = 11, units = "cm", dpi = 600, device = "svg",
         plot = plots.TempTorpor[[m]], scale = 1)
  
  ggsave(filename = file.path(pathWD, paste0(listOutcome[m], "_TempTorpor", ".png")),
         width = 12, height = 11, units = "cm", dpi = 600, device = "png",
         plot = plots.TempTorpor[[m]], scale = 1)
  
  # results table
  tbl.plotResults.TempTorpor <- rbind(tbl.plotResults.TempTorpor, data.frame(rep(listOutcome[m], times = nrow(listResults.TempTorpor.means[[1]])), listResults.TempTorpor.means[[m]]$tempLevel, listResults.TorporTemp.means[[m]]$torporType, listResults.TempTorpor.means[[m]]$torporFlag, listResults.TempTorpor.means[[m]]$emmean, listResults.TempTorpor.means[[m]]$SE, listResults.TempTorpor.means[[m]]$upper.CL - listResults.TempTorpor.means[[m]]$emmean))
}

colnames(tbl.plotResults.TempTorpor) <- c("targetVar", "tempLevel", "torporType", "torporFlag", "mean", "SEM", "CI")
write.csv(tbl.plotResults.TempTorpor, file.path(pathWD, paste0("tblPlotResults", "_TempTorpor", ".csv")), row.names = FALSE)

plots.TempTorpor

```

# save data space image
```{r}
save.image(file.path(pathWD, paste0("RDataSpace_v4", ".RData")))
```

